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Abstract. Multidimensional magneto-hydrodynamical (MHD) simulations coupled with stochastic differential 
■ equations (SDEs) adapted to test particle acceleration and transport in complex astrophysical flows are pre- 

sented. The numerical scheme allows the investigation of shock acceleration, adiabatic and radiative losses as 
well as diffusive spatial transport in various diffusion regimes. The applicability of SDEs to astrophysics is first 
discussed in regards to the different regimes and the MHD code spatial resolution. The procedure is then ap- 
plied to 2.5D MHD-SDE simulations of kilo-parsec scale extragalactic jets. The ability of SDE to reproduce 
analytical solutions of the diffusion-convection equation for electrons is tested through the incorporation of an 
increasing number of effects: shock acceleration, spatially dependent diffusion coefficients and synchrotron losses. 
The SDEs prove to be efficient in various shock configuration occurring in the inner jet during the development 
of the Kelvin-Helmholtz instability. The particle acceleration in snapshots of strong single and multiple shock 



in 

CO 

' acceleration including realistic spatial transport is treated. In chaotic magnetic diffusion regime, turbulence levels 

rfr =< SB 2 > I (B 2 + < SB 2 >) around 0.2 — 0.3 are found to be the most efficient to enable particles to reach the 
highest energies. The spectrum, extending from 100 MeV to few TeV (or even 100 TeV for fast flows), does not 
exhibit a power-law shape due to transverse momentum dependent escapes. Out of this range, the confinement is 
not so efficient and the spectrum cut-off above few hundreds of GeV, questioning the Chandra observations of X- 
ray knots as being synchrotron radiation. The extension to full time dependent simulations to X-ray extragalactic 
?—j | jets is discussed. 
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1. Introduction the physics ijSambruna et alJ l2002|) . The kpc jets show 

nonthermal radio and optical spectra usually associated 
Extragalactic jets in radio-loud active galactic nuclei with sync h ro tron radiation produced by highly relativis- 
(AGN) show distinct, scale dependent structures. At par- tic TeV electrons (positrons may also contribute to the 
sec (pc) scales from the core, superluminal motions have flux )_ The origin of the x _ ray emission is more contro- 
been detected using VLBI technics. The jets decelerate vers i a l and could result from synchrotron radiation or 
while reaching kiloparsec (kpc) scales and power large Inverse Compton (IC) re-processing of low energy pho- 
scale luminous radio lobes. The inner physical condi- tons coming from different sources as synchrotron ra- 
tions are still widely debated. Main uncertainties concern diation (synchro-Compton effect) or cosmic micro-wave 
bulk velocities, matter content, emission and acceleration backgro u nd radiation (CMBR): s ee iMeisenheimer et al 



mechanisms, the way energy is shared between magnetic jS^ , iTavecchio et all lEfflfl . lHarris fc Krawczvnsk: 

field and plasma and finally effects of the turbulent flow ||2002|) for recent reviews 

on relativistic particles. Different acceleration mechanisms have been invoked 

Recent X-ray high resolution observations by Chandra, so far to produce energetic particles, e.g. diffusive 

combined with Hubble space telescope (HST) and radio shock acce l e ration (DSA), second order Fermi accelera- 

data allow unprecedented multi-wavelength mapping of tion in a magnetohydrodynam i c wave turbulenc e (FII) 

the jet structures which lead to improved constraints of fejermann fc StrittmatteJI 19871 iHenri et al.lfl99fl. shock 

drift acceleration (SDA) ijBeeelman fc Kirkl Il990i and 
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magnetic r e conne ction (see for instance iBlackmanl l|l996^ . 
Birk et all l)200l|) and references therein) 1 . With some 
assumptions, all these mechanisms are able to accel- 
erate electrons up to TeV energies and are proba- 
bly at work together in extragalactic jets. Their com- 
bined ^ffectsJmve_OTih^^ 

ever Campeanu fc Schlickeiserl l|l992(h iManolakou et al.1 
(1999)), coupled shock acceleration and spatial trans- 
po rt effects have be en success fully applied to hot spots 
bv lKardashel l|l962h (see also IManolakou fc Kirkl (|2002^ 
and reference therein). Nevertheless, the resolution of 
the full convection-diffusion equation governing the dy- 
namical evolution of the particle distribution function 
does not usually lead to analytical solutions. The parti- 
cle transport and acceleration are closely connected to 
the local magneto-fluid properties of the flow (fluid ve- 
locity fields, electro-magnetic fields, turbulence). Recent 
progress in computational modeling have associated mul- 
tidimensional fluid approaches (hydrodynamical (HD) or 
magnetohydrod ynamical co d es (MHD)) with kine t ic par - 
ticle schemes ll.Tones et al ] I.Tones et al J jlflflfljl . 
iMicono et all l)l999j) h These codes are able to describe the 
effects of shock and stochastic acceleration, adiabatic and 
radiative losses and the results are used to produce syn- 
thetic radio, optical and X-ray maps. The particle trans- 
port is due to advection by the mean stream and turbulent 
flows. In the "Jones et al" approach the shock acceleration 
process is treated using Bohm prescription, i.e. the particle 
mean free path equals to the Larmor radius). These above- 
mentionned treatments neglect spatial turbulent transport 
and introduce spurious effects in the acceleration mecha- 
nism. This leads to an overestimate of the particle accel- 
eration efficiency in jets. 

In this work, we present a new method coupling kinetic 
theory and MHD simulations in multi-dimensional turbu- 
lent flows. We applied the method to the extragalactic 
non- relativistic or mildly relativistic (with a bulk Lorentz 
factor Tjct < 2) jets. Relativistic motions can however be 
handled in case of non-relativistic shocks moving in a rel- 
ativistic jet flow pattern. 

The paper covers from discussions about turbulent trans- 
port and the coupling of kinetic schemes-MHD code to 
more specific problems linked to jet physics. In Section [21 
we review the most important results concerning weak tur- 
bulence theory and exposes the effect of chaotic magnetic 
effects on the relativistic particle (RPs) transport. Section 
El presents the system of stochastic differential equations 
(SDEs) used to solved the diffusion-convection equation 
of RPs. We examine the limits of the SDEs as regards to 
different diffusion regimes and discuss their applicability 
to astrophysics. Section 01 tests the ability of SDEs to de- 
scribe accurately transport and acceleration of RPs in 2D 
versus known analytical results. The MHD simulations of 
jets are presented at this stage to investigate the problem 
of shock acceleration. Section treats RPs transport and 



acceleration in complex flows configurations as found in 
extragalactic jets. We consider the problem of curved and 
non constant compression ratio shocks. We derived analyt- 
ical estimates on the expected particle maximum energy 
fixed by radiative losses or transversal escapes due chaotic 
magnetic diffusivity and MHD turbulence. We report our 
first results on X-ray jets using MHD-SDE snapshots mix- 
ing spatial transport, synchrotron losses, strong single and 
multiple shock acceleration. We conclude in Sectional 

2. Acceleration and spatial transport 

The accurate knowledge of transport coefficients is a key 
point to probe the efficiency of the Fermi acceleration 
mechanisms as well as the spatial transport of RPs in tur- 
bulent sources. We assume a pre-existing turbulent spec- 
trum of plasma waves, retaining the Alfven waves efficient 
to scatter off and accelerate charged particles. The parti- 
cle trajectories are random walks in space and energy, su- 
perimposed to the advection motion induced by the back- 
ground flow, provided that the diffusion time is larger than 
the coherence time of the pitch angle cosine. If the turbu- 
lence level, defined as the ratio of chaotic magnetic com- 
ponents to total one r\ T — < SB 2 > /(< B 2 + SB 2 >) 
is much smaller than unity, the spatial transport parallel 
to the mean magnetic field can be described by the quasi- 
linear theory Before discussing any acceleration mecha- 
nism we shall recall the main results of this theory and 
some of its non-linear developments. 

2.1. Particle transport theories 

During its random walk on a timescale At the position 
of the particle is changed by an amount Axu along 
the mean ordered magnetic field and by Ax± in the 
transverse direction. The ensemble average of both 
quantities vanishes, but the mean quadratic deviations 
are non zero and define the parallel diffusion coefficient 



D|l =< As 



> /2At and the perpendicular diffusion 
coefficient D± =< Ax 2 ^ > /2At. 

For a power-law turbulent spectrum S(k) oc 
Vt (k \max)~^ completely defined by its turbulent 
level tjt, spectral index (3 and maximum turbulent 
scale Xmax, the quasi -linear scatt ering frequency 
v s =< Acos 2 (6>) > /At is l|Jokipiilll96(| 

v, = m n s \fif- 1 pP- 1 . (i) 

Q s is the synchrotron gyro- frequency ZeB /jm^c for a par- 
ticle of charge Ze, mass to* and Lorentz factor 7 and pitch- 
angle cosine p, = cos 9. The Larmor radius Tl = v/£l s and 
the particle rigidity p — lixr^j X m ax- 

The scattering time r s is the coherence time of the pitch- 
angle cosine and can be related to the pitch-angle fre- 
quency v s by t s ~ 1/zAj since the deflection of the pitch- 
angle typically occurs on one scattering time. The quasi- 
linear diffusion coefficients are 



1 lOstrowskH teOOtt) and references therein considered the ef- 
fect of tangential discontinuities in relativistic jets. 
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Nevertheless, this simple approach does not take into ac- 
count the displacement of the guide-centers of particle 
trajectories. When magnetic turbulence is occurring, the 
magnetic field lines are also diffusing, which will amplify 
the transverse di ffusion of par ticles following these mag- 
netic field lines l)Jokit)iil [l969). Indeed including this ef- 
fect in the diffusion dynamics leads to a new transverse 
diffu sion regime, namely the chaotic transverse dif fusion 
iRechester fc R.osenhluthl|T973 iRax fc WhitellTflfll. The 
work done bv lCasse et all ()2002|) presents extensive Monte- 
Carlo simulations of charged particles in a magnetic field 
composed of a regular and a turbulent part, calculated as- 
suming power-law spectra of index (3 (as in Kolmogorov or 
Kraichnan theories). The authors present, using averaged 
spatial displacements over time intervals, the behavior of 
the spatial diffusion coefficients as a function of the parti- 
cles energies as well as turbulence level t\t- The diffusion 
coefficient along the mean magnetic field displays ener- 
getic dependence similar to the quasi-linear theory but 
on any turbulence level. On the other hand, the diffu- 
sion coefficient transverse to the mean magnetic field is 
clearly in disagreement with neo-classical prediction (see 
Eq. [2J • The chaotic transverse diffusion regime is occur- 
ring when the turbulence level is large but can probably 
be extended to lower turbu l ent le vels , as first imagined 
bv IRechester fc Rosenbluthl l)l978)h In lCasse et all l|2002|) 
this regime was observed for all turbulence levels down to 
t\t = 0.03. The resulting transverse coefficient is reduced 
to D± cx Dm with a proportionality factor only depending 
on the turbulence level, namely 



.Dm cx p ' , 

D± cx r)T 3 c\ max p 2 ~ !i 



(3) 



In this paper we will use the above prescription as, unless 
very low t]t, the chaotic diffusion always dominates. 

2.2. Acceleration processes 

In a diffusive shock 2 particles able to resonate with wave 
turbulence, undertake a pitch-angle scattering back and 
forth across the shock front gaining energy. The finite 
extension of the diffusive zone implies some escapes in 
the downstream flow. The stationary solution for a non- 
relativistic shock can be written as f(p) cx p~( 3 + T "^/ T '^\ 
In a strong shock the acceleration timescale T nr . r . ex- 
actly balances the particle escape time scale r esc ijDrurvl 
The acceleration timescale, for a parallel shock is 
TaccDSA = 3/(r — 1) t r , where r = u u /ud is the shock com- 
pression ratio (u u and Ud are respectively upstream and 

2 The shock drift acceleration mechanism has been ap- 
plicd to electron acce l eratio n in extragalactic radio sources by 



lAnastiadis fc Vlahosl il 9931) and references therein. This effect 
will not be considered in the simulations and is not further 
discussed 



downstream velocities of the fluid in the shock frame) and 
t r = (c/ud) 2 r s is the downstream particle residence time. 
The MHD turbulence, especially the Alfven turbulence, 
mainly provokes a diffusion of the particle pitch-angle. 
But the weak electric field of the waves 8E/5B = V a /c 
also accelerates particles. The momentum diffusion is of 
second order in terms of Fokker-Planck description and 
the acceleration timescale is T acc pu — (c/V a ) 2 t s . Note 
that even if the stochastic acceleration is a second order 
process, T aC cFii may be of the same order as T acc SDA in 
low (sub-alfvenic) velocity flows or high Alfven speed me- 
dia as remarked b vlHenri et alJ |l999jh 

In radio jets (see iFerraril l |l985|) arid iFerraril l|l998l) for 
reviews of jet properties) equiparition between magnetic 
fields and non-thermal, thermal plasmas lead to typical 
magnetic fields B ~ I0~ 5 / -4 Gauss, thermal proton den- 
sity n p ~ I0~ 2// ~ 5 cm -3 and thus to Alfven speeds V a /c 
between 7 x I0~ 4 - 0.2. In light and magnetized jets, the 
second order Fermi process can be faster than diffusive 
shock acceleration. We decided to postpone the effect of 
second order Fermi acceleration to a future work. In this 
first step, we mostly aim to disentangle the diffusive shock 
acceleration process, the turbulent spatial transport and 
radiative losses effects shaping the particle distribution. 
We will therefore only consider super- Alfvcnic flows here- 
after. 

3. Numerical framework 

In this section, we present the multidimensional stochastic 
differential equations system equivalent to the diffusion- 
convection equation of RPs. 

3.1. Stochastic differential equations 

The SDEs are an equivalent formulation of the Fokker- 
Planck equations describing the evolution of the distribu- 
tio n fu n ction of a particle population. It has been shown 
bvlltcl <195ll) that the distribution function / obeying 
Fokker-Planck equation as 



21 

at 



" d 



^ — ( Al (t,X)f(t,X)) 

i—1 1 



J2B ik (t,X)Bl(t,X)f(t,X)) (4) 



x N N g2 I \ 

2 IJX/ dX i dX j , 

i=l j=l 3 \k=l / 

at a point X of phase space of dimens ion N, can also be de- 
scribed as a set of SDEs of the form feriills fc Achterberd 
1994) 



N 



d ^ -Mt,x t ) + Y^ BlJ {t,x t ) dWt <> 



dt 



3 = 1 



dt 



i= 1, ..,JV(5) 



3 Ivan der Swaluw fc Achterberg I i200 lT) have investigated 
the coupling between 2D Hydrodynamical code and SDEs 
adapted to the nonthermal X-ray emission from supernova 
remnants. 
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where the Wtj are Wiener processes satisfying < W >= 
W and < (W - W ) 2 >= t - t a (W is the value of 
W at t a ). The diffusion process described by Fokker- 
Planck equations can be similarly taken into account if 
dWi/dt = £i is a random variable with a Gaussian condi- 
tional probability such as 



P(*,el*o,0 



1 



V2ir(t - t a ) 



exp 



2{t-t ) 



(G) 



The Fo kker-Planck e quation governing this population 
will be l|Skillinglll97l 



§j- = -(u-V^+^V-u^ + V^CAiV,/) 
1 9 f d f \ 



(7) 



where Dy is the spatial diffusion tensor and D pp de- 
scribes energy diffusion in momentum space. The term 
stands for synchrotron losses of the electrons. Its ex- 



^syn 

pression is 



L syn 



<7 T B 2 



(8) 



where <tt is the Thomson cross-section. This term can 
easily be modified to account for Inverse Compton losses. 
In term of the variable F — Rp 2 f, these equations can 
be written in cylindrical symmetry (R varies along the jet 
radius and Z along the axial direction) 



dF 
~dt 
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(9) 



Note that this rewriting of the Fokker-Planck equation is 
valid only if Rp > 0. Assuming that the diffusion tensor is 
diagonal, it is straightforward to get the SDEs coefficients. 
These equations can then be written as 



dR 

~dt 
dZ_ 

dt 



= U R 



U, 



dD 



BR 



D 



RR 



OR 
dDzz 



R 

dW z 



dp p_ 
— = — - V • u 

dt 3 



dZ dt 
1 dp 2 D, 



dW R 
' dt 

^JlDzz 



RR 



P 



dp 



Q>synP 



+ ^f^/2Dp~p 



(10) 
(11) 

(12) 



where Ur/z stand for the radial and axial component of 
fluid velocity field. The W are stochastic variables de- 
scribed previously. They are computed using a Monte- 
Carlo subroutine giving a random value £ with zero mean 
and unit variance so that we can build the trajectory of 



one particle in phase spac e from time tk to tk+i = ifc + A£ 
I Marcowith fc Kirklll999^ 
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1 dRD R R 
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(13) 



At + ZzV^DzzAt . (14) 



1 dp 2 Dp 



dp 



a syn p 2 ) At 



£,p^/2DppAt . 



(15) 



It is noteworthy that these algorithms derived from SDEs 
are only valid if the particles are not at the exact location 
of the jet axis, otherwise an unphysical singularity would 
occur. The coupling between the SDEs and a macroscopic 
simulations clearly appears here. The macroscopic simula- 
tion of the jet, using magnetohydrodynamics, would give 
the divergence of the flow velocity as well as the strength 
and the orientation of the magnetic field at the location of 
the particle. Indeed, as shown in the last paragraph, the 
spatial diffusion of particles is mainly driven by the micro- 
scopic o ne, namely by the magnetic turbulence. Since the 
work of I Casse et all ^2^02.) , the behavior of the diffusion 
coefficients both along and transverse to the mean mag- 
netic field are better known. They depend on the strength 
of the mean magnetic field, on the particle energy, and on 
the level of the turbulence i]t- Once the diffusion coeffi- 
cients are known, the distribution function is calculated at 
a time t at the shock front by summing the particles cross- 
ing the shock between t and t+At. The distribution func- 
tion can in princinple be calculated everywhere if statistics 
are good enough. We typically used 5 x 10 5 — 10 6 particles 
per run. 

3.2. Constrains on SDE schemes 
3.2.1. Scale ordering 

Particles gain energy in any compression in a flow. A com- 
pression is considered as a shock if it occurs on a scale 
much smaller than the test RPs mean free path. The ac- 
celeration rate of a particle with momentum p through the 
first-order Fermi process is given by the divergence term 
in Eq.JTJJl, e.g. < ^ >= — § V.v. The schemes used i n 
the present work are explicit l|Kriills fc Achterberdll99^) . 
i.e. the divergence is e valuated at the starting position 
x(tk). Implicit schemes l|Marcowith fc KirlJll999j) use the 
velocity field at the final position x(tk+i) to compute the 
divergence as (u(t k+ i) - u(t k ))/(x(t k+1 ) - x(t k )). 
The particle walk can be decomposed into an advective 
and a diffusive step evaluated at t k and incremented to 
the values R(t), Z(t),p(t) to obtain the new valu es at 
tfc+i . As demonstrated b v lSmith fc Gardiner I lll989h and 
iKloeden fc Platen! lll99ll) it is possible to expand the Ito 
schemes into Taylor series to include terms of higher or- 
der in At and in turn to improve the accuracy of the 
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algorithms. However, both because higher order schemes 
need to store more data concerning the fluid (higher order 
derivatives) and the schemes have proved to accurately 
compute the shock problem in ID, we only use explicit 
Euler (first order) schemes in the following simulations. 
Hydrodynamical codes usually smear out shocks over a 
given number of grid cells because of (numerical) viscos- 
ity. The shock thickness in 2D is then a vector whose 
components are AAT s hock = ((av AR),(a z AZ)), where 
(Ai?, AZ) describes one grid cell and the coefficients 
(a r ,a z ) are typically of the order of a few. We can 
construct, using the same algebra, an advective AX a d v 
and a diffusive AXdifr ve c tors s teps from equations (10- 
P7|l. iKrulls fc Achterberel (|l994^ 

have found that a SDE 
scheme can correctly calculate the effects of ID shock ac- 
celeration if the different spatial scales of the problem sat- 
isfy the following inequality 



AX a d v -C X. 



shock 



< AX diS 



(16) 



In 2D this inequality must be fulfilled by each of the vector 
components, e.g. 



Ai? adv < a r AR < Ai? diff 
AZ adv «a 2 AZ< AZdiff 



(17) 



These are the 2D explicit schemes conditions for the com- 
putation of shock acceleration. The two previous inequal- 
ities impose constrains on both the simulation timcscale 
At, and the diffusion coefficients Drr and Dzz- 

3.2.2. Minimum diffusion coefficients 

The finite shock thickness results in a lower limit on the 
diffusion coefficient. The condition AAT a d v <C AX,ag im- 
plies a maximum value for the time step A^sde that can 
be used in the SDE method given a fluid velocity u (omit- 
ting for clarity the terms including the derivatives of the 
diffusion coefficients) 



At n 



AXjhock/M 



(18) 



In 2D we shall take the minimum Ai max thus derived from 

Inserting this time step into the second part of the re- 
striction, AX a dv <C AXdis yields a minimum value for 
the diffusion coefficient: 



D„ 



l\u\ X 



hock 



(19) 



If the diffusion coefficient depends on momentum, this 
condition implies that there is a limit on the range 
of momenta that can be simulated. The fact that the 
hydrodynamics sets a limit on the range of momenta 
may be inconvenient in certain applications. One can 
in principle circum vent th i s pro b lem by u s ing a daptive 
mesh refinement l)Bergerl Il986l iLeveaud Il998j) . This 
method increases the grid-resolution in those regions 
where more resolution is needed, for instance around 
shocks. The method is more appropriate than increasing 



the resolution over the whole grid. 

Another possibility would be to sharpen artificially shock 
fronts or to use an implicit SDE scheme. This approach 
can be useful in one dimension, but fails in 2D or 3D when 
the geometry of shock fronts becomes very complicated, 
for instance due to corrugational instabilities. 



3.2.3. Comparisons with other kinetic schemes 

As already emphasized in the introduction, up to now, few 
works have investigated the coupling of HD or MHD codes 
and kinetic transport schemes (mainly adapted to the jet 
problem). T he simulations per fo rmed by T. Jo n es an d 
co-workers [[.Tones et alJ l)2002l) . iTregillis et all l|200ll) . 
I.Tones et"al~l ( 19991) 1 present 2 and 3-dimensional synthetic 
MHD-kinetic radio jets including diffusive acceleration at 
shocks as well as radiative and adiabatic cooling. They 
represent a great improvement compared to previous sim- 
ulation where the radio emissivity was scaled to local gas 
density. The particle transport is treated solving a time- 
dependent diffusion-convection equation. The authors dis- 
tinguished two different jet regions: the smooth flows re- 
gions between two sharp shock fronts where the leading 
transport process is the convection by the magneto-fluid 
and the shock region where the Fermi first order takes 
over. This method can account for stochastic acceleration 
in energy, but the process have not been included in the 
published works. The previous distinction relies on the as- 
sumption that the electron diffusion length is smaller than 
the dynamical length as it is the case for Bohm diffusion 
(see discussion in section 13.2. It . However, Bohm diffusion 
is a very pe culiar regime appea ring for a restricted rigidity 
ranges fsee ICasse et all l)2002l) and Eq. |3J). The magnetic 
chaos may even completely avoid it. It appears then essen- 
tial to encompass diffusive spatial transport within MHD 
simulations. 

iMicono et alJ l|l999l) computed the spatial and energy 
time transport of Lagrangian cells in turbulent flows gen- 
erated by Kelvin-Helmholtz (KH) instability. The ener- 
getic spectrum of a peculiar cell is the soluti on of a spa- 
tially averaged diffusion-convection equation I Kardashevl 
1962|). This approach accounts for the effect of fluid tur- 
bulence on the particle transport but suffers from the 
low number of Lagrangian cells used to explore the jet 
medium. The SDE method has the advantage to increase 
considerably the statistics and to allow the construction 
of radiative maps. As the particles are embedded in the 
magnetized jet both macroscopic (fluid) and microscopic 
(MHD waves, magnetic field wandering) turbulent trans- 
port are naturally included in the simulations. 



4. Testing coupling between MHD and SDEs 

So far, particle energy spectra produced by SDEs were cal- 
culated usin g one dimensional prescribed velocity p rofiles 
(see however fvan der Swahiw fc Achterberg I l 200lh V The 
prescriptions described plane shocks as a velocity discon- 
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tinuity or as a smooth velocity transition. The aim of this 
section is to present energetic spectra arising from shocks 
generated by macroscopic numerical code. The tests will 
increasingly be more complex including different effects 
entering in particle transport and acceleration in extra- 
galactic jets. 

In the first subsection, we present very elementary tests 
devoted to control the accuracy of the particle transport 
by SDEs in a cylindrical framework. In the second part, 
we move to jet physics and present the MHD jet simu- 
lations and discuss the results of particle acceleration in 
near-plane shocks produced by Kelvin-Hclmholtz instabil- 
ities. 



4.1. Testing 2D cylindrical diffusive transport 

Before proceeding to any simulations where both MHD 
and SDE are coupled, we have tested the realness of 
our description of the spatial transport of relativistic 
particles. Testing SDEs has already been addressed by 
iMarcowith fc Kirkl ljl999|) and references therein but only 
in a one-dimensional framework. They successfully de- 
scribed the particles acceleration by thin shocks as well 
as the synchrotron emission occurring in the case of rel- 
ativistic electrons. Here, we have added a second spatial 
SDE, for the radial transport, where extra-terms appear 
because of the cylindrical symmetry. One way to test the 
2D transport is to compute the confinement time of a par- 
ticle set in the simple case of a uniform jet with uniform 
diffusion coefficient Drr and Dzz- Let assume we have 
a set of N particles at the jet axis at t = 0. The diffu- 
sion process will tend to dilute this population in space 
and after some time, most of the particles will leave the 
plasma column. Indeed, the average position will be the 
initial position but the spatial variance of these particles 
at time t will be \j2Dt. For the specific problem of a cylin- 
drical jet, let consider a cross section in the Cartesian X 
and Y directions while Z is along the jet axis. The set of 
particles will stop to be confined once 



R 2 jet < 2D xx t + 2D YY t 



(20) 



where Rj e t is the jet radius and the diffusion coefficients 
Dxx and Dyy can be related to Drr by 

< AR 2 > 



Drr 



2At 

Dxx = D Y y 



/ (XAX + )'AY) 2 



R 2 



1 

2At 



(21) 



In this relation, X and Y are two uncorrelated variables (< 
AX AY >= 0). It is then easy to see that the confinement 
time of a set of particles inside a jet is 



R %t 
4D RR 



(22) 



when one consider an infinitely long jet (no particle escape 
in the Z direction) . We have performed a series of calcula- 
tions dealing with one million particles injected near the 
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Fig. 1. Plot of the distribution function F = Rf mod- 
elized by SDE in the case of a uniform spatial diffusion, for 
a fixed Z versus the radial coordinate in jet radius unit. 
The solid curve is the analytical solution obtained from 
Fokker-Planck equation Eq. (0) which is in good agreement 
with computations using SDEs. 
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0.0125 
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0.025 
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9.94 


0.075 
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0.15 
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Table 1. Computations of confinement time T c f for dif- 
ferent diffusion coefficient values and theoretical value of 
this confinement time. Note that the agreement is good 
as far as the confinement time is large. Indeed, the time 
step At = 5 x 1CP 3 to compute them is the same for the 
three runs which leads to different ratio T c f/ At. If this ra- 
tio is too small, the time step is not appropriate to nicely 
modelize the particle transport. 



jet axis with different values of the radial diffusion co- 
efficient. We have set a time step of At = 5.ICP 3 and 
integrated the particles trajectories using the numerical 
scheme EQ. (|15fl . When a particle has reached the jet sur- 
face (i? = Rjet) j we stop the integration and note its con- 
finement time. Once all particles have reached the jet sur- 
face, we calculate the average value of the confinement 
time. In Tab. 14.11 we present the result of the different 
computations. The good agreement between the numerical 
and the estimated confinement times is a clue indicating 
that the spatial transport of the particles in the jet is well 
treated as far as the time step is small enough to mimic 
the Brownian motion of particles. Another way to test 
SDEs in this problem is to look at the distribution func- 
tion of these particles since the analytical solution to the 
diffusion with uniform coefficients is known. The Fokker- 
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Planck equation, in the case of a uniform spatial diffusion 
without any energetic gains or losses, is 



dt 



D 



RR 



d 



R OR 



R dR 



Dzz 



dZ 2 



The radial dependence of f arising from this equation is, 
for an initial set of particles located at the jet axis, 



f(R, z, t)« 



1 



exp 



R 2 \ 



(24) 



On Fig. ^ we plot the distribution function F = Rf ob- 
tained for a set of 5 x 10 5 particles located initially very 
close to the jet axis. The plot is done at a given time 
t = 2 and with D RR = D± = 0.1. The symbols represent 
the numerical values obtained using SDEs while the solid 
line represents the analytical solution from Eg. 1)24(1. The 
good agreement between the two curves is a direct con- 
firmation that the transport of particles is well modelizcd 
by SDEs. 

4.2. MHD simulations of extragalactic jets 

In order to describe the evolution of the jet structure, we 
have e mploy ed the Versatile Advection Code (VAC, see 
iTothl l|l996j) and |http :"7? www I ph ys . uu . nl/| ^toth) . We 
solve the set of MHD equations under the assumption of 
a cylindrical symmetry. The initial conditions described 
above are time advanced using the conservative, second 
order ac curate Total Variation Diminishing Lax-Friedrich 
scheme l|T6th fc Odstrcill 996) with minmod limiting ap- 
plied on the primitive variables. We use a dimensionally 
unsplit, explicit predictor-corrector time marching. We en- 
force the divergence of the magnetic field to be zero by 
applying a projection sch eme prior to every time step 
llBrackbill & Barneslll980h . 

4.2.1. MHD equations 

We assume the jet to be described by ideal MHD in an ax- 
isymmetric framework. This assumption of no resistivity 
v m has consequences on the particle acceleration since the 
Ohm law states the electric field as E = u x B. This elec- 
tric field will vanish in the fluid rest frame so that no first- 
order Fermi acceleration can be achieved by E. In the case 
of a resistive plasma, the electric field (E = B x u + i/, Tl J, J 
density current) cannot vanish by a frame transformation 
and a first-order Fermi acceleration will occur. In order to 
capture the dynamics of shocks, the VAC code has been 
designed to solve MHD equations in a conservative form, 
namely to insure conservation of mass, momentum and 
energy. The mass conservation is 



^ + V-(pu p ) = 0, 



(25) 



where p is the density and u p is the poloi'dal component 
of the velocity. The momentum conservation has to deal 



with both thermal pressure gradient and MHD Lorentz 
force, namely 



dpu 
( 23 ) ~~dT 



V • ^upu 



BB 

Po 



V( 



2/i 



P) =0 , 



(26) 



where P stands for thermal plasma pressure. The induc- 
tion equation for the magnetic field is 



- = -V.(uB-Bu), 



(27) 



The last equation deals with the energy conservation. The 
total energy 



2 



B_ 

2p Q 



P 



r - i 



where T = C P /C V 
governed by 



(28) 

5/3 is the specific heat ratio, is 



ue 



BB 

Po 



U + U 



P 



B^ 

2p Q 



= 



(29) 



In order to close the system of MHD equations, we assume 
the plasma as perfect gas. Thermal pressure is then related 
to mass density and temperature as 



P 



pT 



(30) 



where 5ft is the perfect gas constant and /i p the plasma 
mean molecular weight. 

By definition, these simulations are not able to describe 
microscopic turbulence since MHD is a description of the 
phenomena occurring in a magnetized plasma over large 
distance (typically larger than the Debye distance to in- 
sure electric charge quasi- neutrality) . So, in the case of dif- 
fusion coefficients involving magnetic turbulence, wc shall 
have to assume the turbulence level r]x- 

4.2.2. Initial conditions and boundaries 

We consider an initial configuration of the structure such 
as the jet is a plasma column confined by magnetic field 
and with an axial flow. We add to this equilibrium a radial 
velocity perturbation that will destabilize the flow to cre- 
ate Kelvin-Helmholtz instabilities. The radial balance of 
the jet is provided by the opposite actions of the thermal 
pressure and magnetic force 



B Z (R, Z,t = Q) 
Br(R, Z,t = 0) 

Bg(R, Z,t = Q) 
P(R, Z,t = Q) 



1 , 
. 



(R/Rc) 
1 + (R/R c ) 2 ' 
1 

(l + (R/RcW 



(31) 



where R c is a parameter controlling the location of the 
maximum of Bg and /3 P = 2/i P / B 2 is the ratio of thermal 
to magnetic pressure at the jet axis. All magnetic field 
components are here expressed in B Q units (see next for a 
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definition) . 

The FRI jets are partially collimated flows where some 
instabilities seem to perturb the structure of the jet. Thus 
we will assume in our simulation that the thermal pressure 
is not negligeable in the jet and that the flow is prone to 
axisymmetric Kelvin-Helmholtz instabilities. Thus we will 
assume values of (3 P larger than unity. The sonic Mach 
number is implemented as 



u z (R, Z,t = 0) 

Co 



= M.= 



cosh((i?/i? ) 8 ) 



(32) 



where C is the sound speed at the jet axis. This sound 
speed can be related to jet velocity Ujet by the parameter 
M a = Uj e t/C . This parameter is chosen to be larger than 
one, as the jet is expected to be super-fastmagnetosonic. 
In our simulations, we have M Q = 10 and ug = 0. The 
perturbation that can provoke KH instabilities must have 
a velocity component perpendicular to the flow with a 
wave vector parallel to the flow (e.g. iBodo et all l|l994|) '). 
We have then considered a radial velocity perturbation as 



u R (R, Z,t = 0) = SM J2k'=i sm(kZ2n/L ) 



Co 



exp (5(Z - Ro) 2 ) 



(33) 



where SM a is smaller than unity in order to create a sub- 
sonic perturbation and L is the vertical length of the box. 
The density of the plasma is set as 



p(R, Z,t = 0) 



0.8 



0.2 



(34) 



^cosh((i?/i? ) 8 ) 
where p = T/3 p B 2 /2fi C 2 is the density at the jet axis. 

Physical quantities normalization: 

Lengths are normalized to the jet radius R at the ini- 
tial stage. The magnetic field physical value is given by 
B a while velocities are scaled using sound speed C a = 
Uj et /Mo intimately related the observed jet velocity. Once 
physical values are assigned to the above-mentionned 
quantities, it is straightforward to obtain all the other 
ones. The dynamical timescale of the structure is ex- 
pressed as 

To = ^ = 3-25 x 10 V Mo ( ( Ujet , 1 1 (35) 
° C y \100pcJ \3000km/sJ y ' 

Note that for the single MHD simulations, the evolution 
of the structure does not depend on these physical 
quantities but only on the parameters (3 P , M a , R c , 5M . 
Nevertheless, for each MHD-SDE computation, the 
physical values are injected into SDE equations to derive 
the velocity divergence and the synchrotron losses. 

The MHD simulation are performed using rectangular 
mesh of size 104 x 204 cells, with two cells on each side 
devoted to boundary conditions. The left side of the grid 
(R = 0) is treated as the jet axis, namely assuming sym- 
metric or antisymmetric boundaries conditions for the set 
of quantities (density, momentum, magnetic field and in- 
ternal energy). The right side of the box is at R = AR Q 



and is consistent with free boundary: a zero gradient is set 
for all quantities. For the bottom and upper boundaries 
(respectively at Z = and Z = 8R a , we prescribe peri- 
odic conditions for all quantities, so that when a particle 
reaches one of these regions, it can be re-injected from 
the opposite region without facing artificial discontinuous 
physical quantities. 



4.2.3. Inner-jet shock evolution 

Kelvin-Helmholtz instability is believed to be one of the 
source of flow perturbation in astrophysical jets. The evo- 
lution of this mechanism has been widely i nvestigated ei- 
ther i n a hydrodynamical framework (e.g. iMicono et alJ 
( 200CJ) and refere nce therein) or mo r e rece ntly using MHD 
framework (see iBatv fc Keppensl l|2002|) and reference 
therein). The growth and formation of shock as well as 
vortices in the jet core depend on the nature of the jet 
magnetized or not) and on the magnetic field strength 
Malagoli et aHooA iFrauk et al Jll99fil ljo7es et rtfwM 



IKeppens fc T6thlll999j) . In the particular case of axisym- 
metric jets, it has been shown that the presence of a weak 
magnetic field significantly modifies the evolution of the 
inner structures of vortices. 

We present on Fig. El the temporal evolution of a typi- 
cal inner-jet shock obtained from our computations. After 
the linear growth of the instability (up to t — 19t ), 
the structure exhibits a curved front shock inclined with 
respect to the jet axis. In the frame of the shock, the 
flow is upstream super-fastmagnetosonic, and downstream 
sub-fastmagnetosonic. On both side of the shock, the 
plasma flow remains superalfvenic. This shock configu- 
ration is consistent with a super-fast shock. Rankine- 
Hugoniot relations show that, at a fast-shock front, the 
magnetic field component parallel to the shock front is 
large r in the downstream medium than in the upstream 
one l|Fraix-Burnet fc Pelletierlll99ll) . In the present ax- 
isymmetric simulations, the bending of the poloi'dal mag- 
netic field lines occurring at the shock front creates a lo- 
cally strong Lorentz force that tends to push the structure 
out of the jet. As seen on the following snapshots of FigGJ 
the shock front rapidly evolves toward a plane shape. This 
quasi-plane shock structure remains stable for several time 
units before being diluted. 

4.2.4. Macroscopic quantities 

The SDEs coupled with the MHD code provide approxi- 
mate solutions of the Fokker-Planck equation using macro- 
scopic quantities calculated by the MHD code. Indeed, 
flow velocity and magnetic field enter the kinetic trans- 
port equation and there is no way to treat realistic case 
in astrophysical environments but to model them from 
macroscopic multi-dimensional simulations. Nevertheless, 
one difficulty remains since MHD (or HD) simulations 
only give these macroscopic quantities values at discrete 
location, namely at each cells composing the numerical 
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Fig. 2. Temporal evolution (in r G unit) of a typical inter- 
nal shock occurring within the jet. The grey-scales repre- 
sent density levels (dark for low density and white for high 
density) while solid lines stand for polo'idal magnetic field 
lines. The parameters of this simulations are (3 P = 10, 
M a = 10, R c — 1 and SM = 0.1. This shock arises 
from a initial setup prone to axisymmetric MHD Kelvin- 
Helmholtz instabilities. In the early stage of the shock 
evolution, the shock front displays a bow-shock shape but 
as the simulation goes on, the shock front evolves toward 
a front shock. 



mesh. Hence, these values are interpolated from the grid 
everywhere in the computational domain. If the domain 
we are considering is well-resolved (large number of cells 
in each direction), a simple tri-linear interpolation is suf- 
ficient to capture the local variation of macroscopic quan- 
tities. When shocks are occurring, the sharp transition in 
velocity amplitude is more difficult to evaluate because 
shocks are typically only described by few cells. Thus, the 
calculus of velocity divergence must be done accurately. 
We adopt the following procedure to calculate it: shocks 
are characterized by very negative divergence so at each 
cells we look for the most negative result from three 
methods 



V-u(i,j) 



uz(i,j) -u z (i,j± 1) 



\Z(j)-Z(j±l) 
uz(i,j + 1) - uz[i, j - 1)" 



Z(j + 1) - Z(j - 1) 

R(i)u R {i,j) - R{i ± l)u R {i ± 1, j) 



R(i)(R(i) -R(i±l)) 
R{i + l)u R (i + 1, j) - R(i - l)u R {i - 1, j) 
R(i)(R(i + l)-R(i-l,j)) 



(36) 



This approach ensures that the sharp velocity variation 
occurring within a shock is well described and that no 
artificial smoothing is created in the extrapolation of flow 
velocity divergence. At last, note that the location of the 
most negative Vu corresponds to the shock location. The 
measurement of spectra at shock front will then be done 
by looking at particles characteristics passing through this 
location. 

4.3. Realistic plane shock 

In this subsection we address the issue of the production 
of energetic spectra by plane shocks arising from MHD 
simulations. This issue is a crucial test for the relevance 
of SDEs using the velocity divergence defined in Ea. H36|) . 
We stress that all simulations performed in this paper are 
done using test-particle approximation, i.e. no retroactive 
effects of the accelerated particles on the flow are taken 
into account. 

4.3.1. Strong shock energetic spectrum 

We have performed a series of MHD simulations of cylin- 
drical jets subject to Kelvin-Helmholtz instabilities (cf. 
Sect. I4.2|) . We selected the case of a plane shock (quite 
common in the KH instability simulations) propagating 
along the jet with a radial extension up to the jet radius 
(see Fig. [3J- Its compression ratio is r ~ 4 (measured 
by density contrast) and constant along the shock front. 
We have chosen a particular snapshot of the structure dis- 
played on Fig. |3J By rescaling the vertical velocity in order 
to be in the shock frame (where the down and up-stream 
velocities are linked by Wdown — Uup/r), we first consider 
this shock with infinite vertical boundaries and reflective 
radial boundaries. Namely, we set that if the particle is 
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Fig. 3. Zoom in a jet snapshot where Kelvin-Helmholtz 
instabilities are active. The parameters of the MHD simu- 
lations are the same as in FigEl The grey levels represent 
the density levels while the white lines are magnetic sur- 
faces. A shock arises in the core of the jet (R < 1) with 
a plane shape perpendicular to the jet axis. Using a large 
number of particles like the one which trajectory is dis- 
played with a thick white line, we measure, in the shock 
frame, the stationary energetic spectrum of particles at 
the shock front. 



escaping the domain at Z < Z min = or Z > Z max = 8, 
we take the velocity to be u p (Z > Z max ) — u p (R, Z max ) 
(same thing for Z < Z m i n ). The condition allows for parti- 
cles far from the shock to eventually return and participate 
to the shaping of F(p). The reflective radial boundaries are 
located at the jet axis R = (to avoid the particle to reach 
R = where SDEs are not valid) and R — 1. Such bound- 
aries ensure that no particle can radially escape from the 
jet during the computation. The constant value of the dif- 
fusion coefficients Dzz and Dan must fulfill relations (|17|l 
and l|19|l . Actually, in the particular case of a plane shock 
propagating along the vertical axis, only Dzz must ful- 
fill previous relations, namely Dzz > D min — X s j l \uz\/'2- 
The shocks width X s h is defined as the location of the 
most negative velocity divergence of the flow. Typically, 
this width corresponds to the size of a mesh cell in the 
case of strong shock. We can then safely set Dzz — 0.4 
as we will have Dzz = 10D m i n . The radial diffusion coef- 
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Fig. 4. Energetic spectrum of particle population injected 
at momentum p a in the MHD jet of Fig. [3J The spectrum 
is measured at the shock front and is in a very good agree- 
ment with DSA theory predicting a power-law of index —4 
for a single plane shock of compression ratio r = 4. Note 
that in this computation the diffusion coefficients have 
constant values fulfilling relation (|17fl . 



ficient is tuned as Drr = 0.01 and will enable particle to 
explore the shock front structure. On Fig. 01 we display the 
results of the use of SDEs on a particle population injected 
at momentum p = p a and propagating in snapshot repre- 
sented by Fig. |21 We easily see that the resulting spectrum 
is a power-law of index —4 completely in agreement with 
DSA theory (see Section^. The existence of a few particle 
with p < p arises from the fact that outside the shock, 
the velocity divergence is not equal to zero, as it would 
be with a prescribed velocity p rofile ijKrulls fc Achterberd 
I1994 iMarcowith fc Kirklll999h . Note that in the absence 
of other energetic mechanism (as second-order Fermi ac- 
celeration or synchrotron losses), the simulation is inde- 
pendent of the physical value of p a as the diffusion coeffi- 
cient is independent of p. 



4.3.2. Single shock with synchrotron losses 

For electrons, the acceleration occurring within shock may 
be balanced at the cut-off by radiati ye losses due to the 
presence of the jet magnetic field. I Webb et alJ l|l984h 
has presented a complete analytical resolution of Fokker- 
Planck transport equation including both first-order Fermi 
acceleration and synchrotron emission. In particular, they 
show that the energetic spectrum exhibits a cut-off at a 
momentum p* depending on spatial diffusion coefficient 
and velocity of the flow. The choice of the injection en- 
ergy p a c of electrons is determined by the lower boundary 
of the inertial range of magnetic turbulence. Indeed, to 
interact with turbulence and to spatially diffuse, electrons 
must have momentum larger than pi — niiVA, where 
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Fig. 5. Energetic spectrum of energetic population in- 
jected at momentum p a c = lOOMeV in the MHD jet 
of Fig. |2| and subject to synchrotron emission. The spec- 
trum is measured at the sho ck front. It matches the so- 
lution of IWebb et alJ l) 19841) , in particular for the cut- 
off momentum p* where synchrotron losses balance shock 
energy gains. For this computation, we have considered 
magnetic field given by the MHD code. The upstream 
velocity U up = 300km/s, and the mean magnetic field 
B a = 100 uG. 
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Fig. 6. Energetic spectrum resulting from acceleration by 
multiple shocks in an extragalactic jet. The power-law 
arising from this simulation matches exactly the result 
given by DSA theory in the case of no particle escape 
from the shock region (t esc 3> t acc ). To achieve this simu- 
lation, we have considered the snapshot in Fig|3|but with 
vertical re-injection of escaping particles. These bound- 
aries mimic the effect of multiple shocks interaction with 
particle during their propagation along the jet. 



is typically the proton mass and Va is the Alfven speed 
l|Lacombellr977h . The energy threshold corresponds to 

et = Pi c ~ 900MeV {^^j (37) 

In our simulations, we assume an Alfven speed V a ~ 
2.2 10 s cm/s ~ c/100 (see the discussion in Section l2~2~|l 
leading to p D c = lOOMe^ > piC. As previously noted, the 
Alfven speed in extragalactic jets can reach appreciable 
fraction of the light speed. An increase of V a leads to an 
increase of the particle injection threshold and a decrease 
of the dynamical momentum range explored. In that case, 
the Fermi second order effect must be included in our SDE 
system (via the diffusive term in momentum in Eq. I15|) . 
Time dependent simulations (in progress) will include the 
associated discussion. 

The result of the simulation including synchrotron losses is 
displayed on Fig. [SJ When assuming a constant magnetic 
field and diffusion coefficien ts, the cut-off energy e* = p*c 
reads as l|Webb et alJll984^ 

mic 3 2tt 9 r — 1 

c* = _ u z 

D zz a T B2 up r(r + l) 

Fig- El displays the spectrum at the shock in case of a the 
magnetic field obtained from the MHD code. The cut-off 
is in good agreement with the resulting spectrum despite 
the numerical simulation is considering a spatially varying 
magnetic field. The Fig. |5] also shows the case of a con- 
stant magnetic field taken as the average of the previous 
one. 



4.3.3. Multiple shocks acceleration 

The presence of multiple shocks increases the efficiency of 
particle acceleration. In multiple shocks, the particles ac- 
celerated at one shock are advected downstream towards 
the next shock. The interaction area is enhanced, so as 
the escaping time. The general expression of the distribu- 
tion function at shocks front, log / oc — (3 + t acc /t esc ) logp 
will then tend to log/ cx —3 logp. This multiple shocks 
acceleration may oc cur in jets where numer ous internal 
shocks are present l|Ferrari fc Melroselll997h . We intend 
to modelize this effect using the same snapshot as in pre- 
vious calculations but changing the nature of the vertical 
boundaries. Indeed, since we are modelizing only a small 
part of the jet (typical length of 800pc), we can assume 
that if a particle is escaping by one of the vertical bound- 
aries, it can be re-injected at the opposite boundary with 
identical energy. The re-injection mimics particle encoun- 
ters with several parts of the jet where shocks are oc- 
curring. Physical quantities are set to same values than 
in paragraph dealing with single plane shock. The result 
of the simulation is displayed on Fig. and on Fig. [3 
when synchrotron losses are considered. On Fig. El the 
spectrum reaches again a power-law shape but with a 
larger index of —3, consistent with previous statements. 
When synchrotron losses are included in SDEs (Fig. EJ, 
we find a similar spectrum than for single shock but with 
some differences. Namely, the curve exhibits a bump be- 
fore the cut-off. This bump can easily be understood since 
the hardening of the spectrum enables particles to reach 
higher energies where synchrotron losses become domi- 
nant. Thus an accumulation of particles near the cut-off 
momentum p* will occur. The bump energy corresponds 
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Fig. 7. Momentum spectrum of electrons injected at 
p a c — lOOMeV in a jet prone to multiple shocks accelera- 
tion. Note the bump occurring because of the synchrotron 
cut-off that tends to accumulate particles at the cut-off 
momentum p* . 

to the equality between radiative loss timescale and mul- 
tiple shock acceleration timescale. The last timescale is 
larger than the timescale required to accelerate a particle 
at one isolated shock be cause of the advection of particles 
from shock to next one l|Marcowith fc Ki'r^ ll999'). 

5. Acceleration at complex shock fronts 

The shock structures are subject to important evolution 
during the development of the KH instability. We now 
investigate the particle distribution function produced at 
these shocks using the SDE formalism. All the shock ac- 
celeration process here is investigated using snapshots of 
the MHD flow. 

5.1. Plane shocks with varying compression ratio 

In astrophysical and particularly jet environments, (weak) 
shocks occurring within magnetized flows in the early 
phase of KH instability (see Fig- El) are non-planar and/or 
with non-constant compression ratio along the shock sur- 
face. We first consider analytical calculations of the parti- 
cle distribution produced in such shocks that extend previ- 
ous works and we complete our estimates using the MHD- 
SDE system. 

5.1.1. Analytical approach 

The theory of DSA explains the energetic spectrum of dif- 
fusive particles crossing plane shock with constant com- 
pression ratio r. Even when the plane structure is relaxed 
ilDrurvl Jl983)) the compression ratio is usually assumed 
as constant along the shock front. In astrophysical jets, 
complex flows arise from the jet physics so that even the 
plane shock assumption is no longer valid implying a non- 
analytical derivation of the particle distribution function. 
Nevertheless, it seems obvious that if the shock front is 



not strongly bent, the particle acceleration process should 
not be strongly modified. 

Let us first quantify this assertion. We calculate the mean 
momentum gained by a particle during one cycle (down- 
stream — > upstream — > downstream) 



< Ap > 
P 



g(r(*)-l)- 



(39) 



we assume that, during this cycle, the particle sees the 
local structure of the shock as a plane (v is the speed of 
the particle), e.g. the spatial scale where the shock bends 
is large compared to the particle diffusive length. 
Here, contrary to the standard DSA theory, the energy 
gain depends on the location of the shock crossing of the 
particle. The probability for a particle to escape from the 
shock during one acceleration cycle is however still given 
by the usual DSA theory, namely r]k = Aud/vk (wfe is the 
speed of the particle during the fcth cycle) . The probability 
that a particle stays within the shock region after n cycle 
Pr n can be linked to the mean momentum gain after n 
cycle as 



In Pr„ 



Inpn/po J2k=i M 1 + Vk(rk - l)/3) 



(40) 



The compression ratio r depends on the number of the 
cycle since in reality, the particle is exploring the front 
shock because of the diffusion occuring along the shock 
front. This expression can be simplified if we assume the 
flow background velocity very small compared to particle 
velocity (% <C 1 for a non-relativistic shock) and that 
particles are ultra-relativistic (r/k = rj = Aud/c). The ex- 
pression then becomes 

111 7V " • " . (41) 



In Pn/Po 



En 
k=l r k 



The sum of the different compression ratios experienced by 
particle population can be approximated using the aver- 
age compression ratio measured along the shock front r rn . 
Indeed, each particle interacting with the shock are prone 
to numerous cycles of acceleration and then the sum re- 
maining in Ea. (|41|) can be expressed as r k — nr m- 
Hence, the energetic spectrum is a power-law but with 
an index controlled by the mean value of the compression 
ratio all over the shock front, namely 



V f oc 



d In Pr„ 

dp n 



oc p 



-3(r m +2)/(r m -l) 



(42) 



In this demonstration, the compression ratio profile itself 
is not involved in the spectrum shape but only its aver- 
age value r m , as long as one can consider the shock to 
be locally plane . Ea. (|42[) generalizes the results provided 
by iDrurvl l)l983j) concerning curved shocks with constant 
compression ratio. If the plane shock assumption is re- 
laxed, numerical simulations are necessary. 
In order to complete this result, we have performed several 
numerical calculations where a mono-energetic population 
of relativistic particles are injected with momentum p be- 
hind an analytical prescription describing a plane shock 
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Fig. 8. Spectrum produced by plane shocks with spatially varying compression ratio r. The left panel represents 
compression ratio profiles along the shock front while the right panel displays the resulting spectra. The calculations 
1, 2 and 3 have different r— profiles but the same r m — 4. The resulting spectra are corresponding to curve 1,2 and 
3 on the right figure. We can see that this curves are almost the same and are very close to a power-law with index 
— (r m + 2)/(r m — 1). Curves 4 and 5 correspond to compression ratio profile with mean values r m equal to 3.5 and 3. 
The corresponding spectra on the right figure are again consistent with power-law with indices controlled by their r rn . 
All particles are injected at p Q c — 100 MeV along the shock. 



with varying compression ratio (the shocks are test ex- 
amples) . The result of this numerical test is displayed on 
Fig.©. 

In this test, we have done three calculations with three 
different compression ratio profiles (curves 1, 2 and 3) but 
with identical average value r m = 4. Setting both vertical 
and radial diffusion, we have obtained the spectra 1, 2 and 
3 displayed on right panel of Fig.JHJ. These three curves 
are almost the same. On two other calculations, we have 
chosen linear profiles with different values of r m (curves 4 
and 5): again a power-law spectrum is found with indices 
consistent with previous analytical statements. This con- 
clusion is correct only if during on cycle the particle mean 
free path along the shock front is small compared to its 
curvature and if during many cycles the particle is able to 
explore the whole shock structure. 

5.1.2. Locally-plane shock 

The previous considerations can be applied to a non- 
planar shock produced in the early stage of the ax- 
isymmetric Kelvin- Hclmholtz instability. The inner shocks 
tend to evolve from curved fronts in the early phases of 
the instability toward plane shocks, perpendicular to the 
jet axis (see Sect 14 .2. 31 and FigEJ). On FigElthe curvature 
radius of the shock is typically of the order of the jet ra- 
dius while the obliquity angle (between the shock front 
and the jet axis) ranges from zero to 27°. For such a low 
obliquity the shocks are subluminal. A more subtle con- 
sequence of the non-constancy of the compression ratio 
is that the electric fields generated along the shock front 
cannot be canceled by any Doppler boost. In other words, 
complex shocks do not have a unique de Hoffman- Teller 
frame. This problem strongly complexifies the particle ac- 
celeration and transport in jets and is postponed to fu- 



ture works especially treating strongly oblique (or even 
perpendicular) shocks. In the present paper, the MHD 
shocks are only weakly oblique and non-relativistic (the 
effects of electric fields on particle acceleration are ne- 
glected). In principle, once the electro-magnetic field is 
known throughout the jet, the systematic electric effects 
on particle trajectories can be implemented in the SDE 
system. 

Keeping the same prescription for diffusion coefficients 
than in previous section (constant diffusion coefficients 
and radial reflective boundaries) we first have to verify the 
quasi-planar condition of the shock. To this aim, we form 
the ratio of the typical diffusion length occurring during 
one cycle along the shock front (L±) and the curvature ra- 
dius. The duration of one acceleration cycle is controlled 
by the residence time at the shock t res = 2Dzz/u^ (as- 
suming that it is composed of n identical cycle) . The num- 
ber of cycle needed for the particle to escape the shock 
is obtained when the escaping probability after n cycles 
is equal to unity, namely X^fe=i Vk = 4nu ( j/c when par- 
ticles are relativistic. The duration of one cycle is thus 
r = 8Dzz /ud,c- The criterion for considering a shock as 
locally plane will be 

-r- <- <«> 

where -D_l,s is the maximal value of the diffusion coeffi- 
cient in the direction parallel to the shock front. With the 
previously prescribed diffusion coefficients the ratio has 
a maximum value equal to 3.4 10~ 2 . This value is small 
compared to unity which means that during one cycle, 
the particle will interact with a zone of the shock where 
the compression ratio is almost constant. On the other 
hand, this ratio is not so small and within a few cycles of 
acceleration particles will explore a significant part of the 
shock front. 
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Fig. 10. Compression ratio profile (left) and energetic spectrum (right) of accelerated particles by the shock displayed 
on FigEl This shock is not a plane shock and does not have constant compression ratio along its shock front. The 
curvature radius of this shock is much larger than the mean free path of a particle so, locally, the shock can be 
considered as a plane shock. The resulting spectra (with only acceleration or with synchrotron losses included) are 
consistent with a plane shock with compression ratio 2.7 which is close to the average value of the compression ratio 
of this shock, namely r m = 2.68. 



Figure (|10fl shows the energetic spectrum produced in 
such curved shock. The result is close to a power-law of 
index / oc p~ A - 7 and when synchrotron losses are taken 
into account, the cut-off energy corresponds to the case 
of a uniform shock with constant compression ratio equal 
to 2.76. The cut-off, given by Ea. (|38|l is close to the value 
obtained on the plot. We postpone to Section 15.2.21 the 
comparison between particle acceleration timescale and 
shock survival timescale in the different phases of the jet 
evolution. We can however anticipate here that for typical 
jet parameters the former is smaller than the latter. This 
validates our results obtained using MHD snapshots. 



5.2. Strong shock acceleration and spatial transport 

We now consider the shock acceleration and spatial trans- 
port in chaotic magnetic field in strong shocks occurring 
in the late phase of the KH instability w here the most 
efficie nt particle acceleration is expected ()Micono et alJ 
The validity of the snapshot approach is tested 
against the survival of the shocks. We investigate the effect 
of radial escape on the particle distribution in the single 
and the multiple shock configuration. 



5.2.1. Maximum energies expected and electron 
transport 

The maximum electron energy is limited by radiative or 
escape losses. In case of synchrotron radiation, the loss 
timescale is ti oss ~ 1.2 10 4 (-B m G)~ 2 Eq\ v yr which com- 
pared to the acceleration timescales presented above leads 
to electron with energies j ma x ~ 10 8 (Uj et /c) 3 ^ 2 , around 



energy are expressed in mGauss and in GeV units respec- 
tively) . 

In a quite general way, the radial and vertical diffusion 
coefficients can be written as 
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(44) 



where B is the total magnetic field amplitude. The con- 
finement time is driven by the radial diffusion coefficient 
Drr which may be expressed in term of t]t as 



rj T 



(45) 



where Dm = D jr\T and a stands for the average value 
of \Br/B\ all over the simulation box. The coefficient D a 
may eventually depends on the particle momentum. It can 
easily be seen that Ea. (|45|l has a minimum value for 



mm 



1.3 (1 



2)1/2 



1/2.3 



(46) 



At ij m i n , particle confinement reaches its maximum (see 
Eq. |22Jl- Typically, radio jets do not display opening an- 
gle larger than a few degrees leading to r]™ m of the order 
of 0.2. Paradoxically, low turbulence levels do not pro- 
vide efficient confinement since largest diffusion motion 
occurs along the magnetic field which have locally radial 
components. Within a timescale ti oss electrons are able to 
explore distances 



1 TeV for Uj et /c — 0.1 (the magnetic field and the particle AR = \J 4 Drr ti 



(47) 
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Fig. 9. Same plot as Fig.©. The MHD simulation pa- 
rameters are the same as in FigOU except for SM a = 0.5. 
This snapshot is selected during the developing phase of 
the shock where the structure is evolving toward a plane 
shock. The inclined (with respect to the jet axis) part of 
the shock front are affected by a strong magnetic bending 
due to magnetic conservation through the shock. 

In the chaotic magnetic regime, Eq.© leads to 

Ai? ~ 15 [ (^-^^V E&* pc. (48) 

We consider a mean turbulence level t]t — i]jp m , and 
assumed a magnetic field B = 100 fiG, and a maxi- 
mum turbulence scale X max ~ Rj C t- For a ~ 2°, we get 
AR(1 TeV) - 9 pc and Ai?(l GeV) - 90 pc. The high 
energy electrons are only able to explore about one tenth 
of the jet radius and can be considered as confined to 
the region where they have been injected. The GeV elec- 
trons can explore larger fraction of the jet and escapes are 
expected to steepen the particle distribution. These are 
averaged results, AR is sensitive to the magnetic field, for 
example if B decreases (increases) by one order of mag- 
nitude AR increases (decreases) by a factor ~ 14. Along 
the jet, particles are advected from one shock to the next 
on timescales AZ s h oc k/Vfiow, where a mean inter-shock 
distance AZ s hock ~ 1 kpc and Vji ow ~ 10~ 2 / _1 c lead 
to T ac ;,j|| ~ 3 10 4 / 5 yr. The high energy part of the elec- 
tron distribution is then produced by one shock and can 
hardly be re-accelerated in a second one downstream. The 



spectrum at these energies strongly depends on the shock 
compression ratio. At lower energies GeV electrons dis- 
tribution can be subject to either transversal escapes or 
multiple shock effects. For both populations, the electrons 
accelerated at inner shocks remain within 1 kpc from their 
injection points, this clearly separates the inner jet to the 
Mach disc and justifies a fortiori our approach simulat- 
ing only the kiloparsec scale jet. It also clearly appears 
that the spatial transport issue addresses to the multi- 
wavelength morphologies of jets. We know make these 
statements more precise using the coupled MHD-SDE sys- 
tem. 

5.2.2. Single shock 

So far, we have presented numerical calculations using re- 
flective radial boundaries (no particle losses) and constant 
diffusion coefficients. In this section, we choose to remove 
step by step these two constraints. Starting from the snap- 
shot of Fig|2| we first remove the outer reflective boundary 
and consider any particle having R > Rj e t as lost. Then 
we adopt diffusion coefficients given by Eq. © since they 
arise from a transport theory consistent with high turbu- 
lence levels rjT and are confirmed numerically. Quasi-linear 
theory docs likely apply at very low turbulence levels im- 
plying high parallel diffusion coefficient and acceleration 
timescales. Expected spectra must then be softer than the 
same spectrum obtained in chaotic regime. 
First, as an illustration of escape effects, we consider the 
case of constant diffusion coefficients, namely Dzz = 1 
and Drji — 2.10" 2 . The resulting spectrum can be seen 
on Fie llll and is consistent with a harder power-law. In 
previous simulations, the escaping time was defined as the 
time needed by the flow to advect a majority of RPs away 
from the shock. Here the effect of the confinement inside 
the jet if lower than the escaping time from the shock 
will be the main source of particle losses. The distribu- 
tion function reads as log / oc — (3 + t acc /ti oss ) \ogp where 
ti oss = min(t con f,t esc ) and will stay as a power-law as 
long as the escaping time is not momentum dependent. 
In our example t loss = t conf = (R jet - R mo ) 2 /^D RR = 
10.125, where Rmj is the average radius of injected par- 
ticles. The resulting spectrum index on Fie. ll II is in good 
agreement with this estimate since the ratio t a cc/t C onf — 
&D / {t con f u 2 d (r — 1)) = 1-26 and the plot representing the 
spectrum done with these constant diffusion coefficients 
has a power-law index as / oc p~ 4 ' 25 . 
Secondly, we discuss the case of Kolmogorov turbu- 
lence and keep tjt free in order to check its influence 
on the transport of particles. The three last plots in 
Fig ll II represent simulations performed without outer re- 
flective boundaries and diffusion coefficients as described 
by Eq^J The simulations account for energy as well as 
spatially [B r and B z are both function of r and z) depen- 
dent transports. Each curve corresponds to a value of the 
turbulence level r\r = 0.05,0.2 and 0.9. In a Kolmogorov 
turbulence D\\ oc D± oc p 1 ^ 3 , t con f oc oc p~ 1 / 3 which 
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Fig. 11. Spectra produced at the shock front displayed on 
Fig|3| without outer reflective boundary, namely with ra- 
dial particle losses. The upper plot represents a spectrum 
done with constant diffusion coefficients (Dzz = 1 and 
Drr = 0.02). The radial losses modify the spectrum by 
increasing the index of the power-law from —4 to —4.25. 
The last three curves are spectra produced by using real- 
istic diffusion coefficients given by Eq. © . The momentum 
dependence of these coefficients modifies the shapes that 
are no longer power-laws (see Sect l5.2.2*|) . The upper plot 
has an arbitrary normalization unrelated to the three last 
curves. 



leads to a confinement time decreasing while increasing 
momentum and a convex spectrum. At a low turbulence 
level, the ratio t aC c/tconf is large, increases with the parti- 
cle momentum and leads to softer spectra with low energy 
cut-off (at few GeV/c). In order to get significant particle 
acceleration and large energy cut-off (beyond 1 TeV) tur- 
bulence levels rjT > 0.1 seem mandatory. The maximum 
confinement is obtained for turbulence level t\t compati- 
ble with Eq. JUJl. 

One important issue to discuss about is the validity of 
our results while considering snapshots produced from the 
MHD code. It appears from Fig.[3]that both weak curved 
and strong plane shocks survive a timcscalc of the order 
of 5 t . The shock acceleration timescale of a particle of 
energy i?GeV may be exp ressed as T„n r ~ 20 Dzz/U^ h for 
a compression ratio of 4 {5rur3[l983), where U s h ~ 10 c s 
is the shock velocity. Using the Eq. and ^] we end up 
with a typical ratio r acc /ro ~ 10~ 2 rfj} Eq^ v B ir ^l z . Our 
snapshot then describes well the shock acceleration (e.g. 
Tucc/tq < 1) up to TeV energies unless the turbulence level 
is very low and the magnetic field much lower than 100/iG. 
The conclusion is the same for curved shocks as the ac- 
celeration timescale is smaller in that case. However, time 
dependent simulations are required to a more exhaustive 
exploration of the jet parameter space and to test the dif- 
ferent turbulence regimes. 



Fig. 12. Same plot as in Fig^Jbut with periodic verti- 
cal boundaries. This setting mimics the effect of multi- 
ple shocks acceleration. As previously, inclusion of radial 
particle losses affects spectra: a softening of the spectral 
index, cut-off energies dependent on t/t ■ 

5.2.3. Multiple shock-in-jet effects 

The radial losses should also affect the transport of parti- 
cles encountering several shocks during their propagation. 
This description is pending to the possibility of multiple 
strong shocks to survive few dynamical times. This issue 
again requires the time coupling between SDE and MHD 
simulations to be treated. 

However, the general statement about the distribution 
function is still valid but at the opposite of previous multi- 
ple shocks acceleration calculations (see Section l4.3.3|) the 
lack of confinement is the only loss term. On Fig 1 121 we 
have performed the same calculations as in the previous 
paragraph except that we impose periodic vertical bound- 
aries where particles escaping the computational domain 
by one of the vertical boundary is re-injected it at the op- 
posite side keeping the same energy. 

We again start with our fiducial case displaying the spec- 
trum obtained from calculations done with constant dif- 
fusion coefficients, i.e. Dzz = 1 and Dan = 2 10~ 2 (the 
upper plot). The power-law index is modified and equals 
to —3.13 instead of —3 as obtained in calculations without 
radial losses. This result is close to the analytic estimate 
since t acc /t con f ~ 0.11 in that case. In the chaotic diffu- 
sion regime the same behavior is observed in the spectra, 
e.g. convex shape, low energy cut-off at low turbulence 
levels. In this diffusion context, multiple shock accelera- 
tion is again most efficient for t/t ~ 0.2 — 0.3 and tends 
to produce hard spectra up to 10-100 GeV for electrons 
without radiative losses. The spectrum cut-off beyond 10 
TeV. 

In Fig El we have included synchrotron losses effects 
in one of the most favorable case (j]t = 0.25) in the 
chaotic regime. The resulting spectrum shows a charac- 
teristic bump below the synchrotron cut-off lying around 
a few ten GeV. This hard spectrum may be inte rmittent in 
jets as already noticed bv lMicono et al . (1999). The spec- 
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Fig. 13. Multiple inner-jet shocks spectrum including syn- 
chrotron losses for t\t = 0.2, B a = 10/zG and Uj e t = 
c/5. The synchrotron cut-off energy lies at a few TeV. 
The most energetic electrons at the bump of this dis- 
tribution would have a peaked synchrotron emission at 
v ~ 1.9 10 16 Hz which corresponds to UV/X-ray emission 
(hv~ 20eV). 

trum and bump may also not exist because of non-linear 
back-reaction of relativistic particles on the shock struc- 
ture (this problem require the inclusion of heavier par- 
ticles in the simulation). Beyond the electrons loss their 
energy before reaching a new shock as discussed in Section 
15.2.11 The magnetic field used is 10 fiG and suggests that 
higher values are apparently not suitable to obtain TeV 
electrons. The synchrotron peaked emission of the most 
energetic electrons of this distributions gives an idea of 
the upper limit of radiative emission achievable by this 
inner-jet shock. In a 10/iG ma gnetic field, these electrons 
radiate UV/X-ray photons as l)Rvbicki fc Lightmanll 979) 



0.29 



3he 2 eB 



20eV (wrev) 2 ''(£g) 



(49) 



The electron population computed here does not go be- 
yond 50 TeV, which then suggests an energy emission up- 
per limit around ~ 0.5keV. The maximum energy scales 
as Uj e t ( see Ea.(|38[ll and can be significantly increased in 
case of fast jets (with speeds up to c/2 the limit of the 
validity of the diffusion approximation). 
In case of inefficient confinement, e.g. t/t different from 
0.2-0.3, this result also suggests that the synchrotron 
model may in principle not account for the X-ray emis- 
sion of extragalactic jets probably dominated by another 
radiative mechanism (for instance the Inverse Compton 
effect). However again, we cannot draw any firm conclu- 
sions about this important issue and postpone it to the 
next work treating full time dependent simulations. 

6. Concluding remarks and outlook 

In the present work, we performed 2.5D MHD simulations 
of periodic parts of extragalactic jets prone to KH insta- 



bilities coupled to a kinetic scheme including shock accel- 
eration, adiabatic and synchrotron losses as well as appro- 
priate spatial transport effects. The particle distribution 
function dynamics is described using stochastic differen- 
tial equations that allow to account for various diffusion 
regimes. 

We demonstrate the ability of the SDEs to treat multi- 
dimensional astrophysical problems. We pointed out the 
limits (D m i n defined in Eq. I19|l imposed by the spatial 
resolution of the shock on the diffusion coefficient. The 
SDEs are applicable to particular astrophysical problem 
provided D > D m i n . We perform different tests in 2D 
showing consistent results between numerical simulations 
and analytical solutions of the diffusion-convection equa- 
tion. Finally we demonstrate the ability of the MHD-SDE 
system to correctly describe the shock acceleration pro- 
cess during the evolution of the KH instability. Complex 
curved shock fronts with non constant diffusion coeffi- 
cients that occur at early stage of the instability behave 
like plane shock provided the diffusion length is smaller 
than the shock curvature. The equivalent plane shock has 
a compression ratio equals to the mean compression of the 
curved shock. In case of strong plane shocks which develop 
at later stages of the KH instability, we found that the in- 
clusion of realistic turbulent effects, e.g. chaotic magnetic 
diffusion lead to complex spectra. The resulting particle 
distributions are no more power-laws but rather exhibit 
convex shapes linked to the nature of the turbulence. In 
this turbulent regime, the most efficient acceleration oc- 
curs at relatively high turbulence levels of the order of 
~ 0.2 — 0.3. The electron maximum energies with syn- 
chrotron losses may go beyond 10 TeV for fiducial mag- 
netic field values in radio jets of ~ 10/iG and the spec- 
trum may be hard at GcV energies due to multiple shock 
effects. 

However, in this work, SDEs were used on snapshots of 
MHD simulations neglecting dynamical coupling effects, 
preventing from any complete description of particle ac- 
celeration in radio jets. Such dynamical effects encompass 
temporal evolution of shock, magnetic field properties and 
particle distribution. The time dependent simulations will 
permit us to explore the parameter space of the turbu- 
lence and to critically test its different regimes. 
The simulations have also been performed in test-particle 
approximation and do not account for the pressure in RPs 
that may modify the shock structure and the acceleration 
efficiency. This problem will be addressed in a particular 
investigation of shock-in-jet acceleration including heav- 
ier (protons and ions) particles. Nevertheless the present 
work brings strong hints about the ability of first order 
Fermi process to provide energetic particles along the jet. 
Our first results tend to show that synchrotron losses may 
prevent any electron to be accelerated at high energies re- 
quiring cither supplementary acceleration mechanisms or 
other radiative emission processes to explain X-ray emis- 
sion as it has been recently claimed. Future works (in 
progress) will account for these different possibilities. 
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